package com.qyx.pmpucat.util;// 导入所需的 Apache Commons Math 库中的类

import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.integration.TrapezoidIntegrator;
import org.apache.commons.math3.distribution.NormalDistribution;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.RealMatrix;

public class Grm {

    public static final double D = 1.0;

    // 定义一个内部类 PiResult，用于存储 Pi 及其导数的结果
    public static class PiResult {
        public double[][] Pi;  // 概率矩阵
        public double[][] dPi;  // 一阶导数矩阵
        public double[][] d2Pi;  // 二阶导数矩阵
        public double[][] d3Pi;  // 三阶导数矩阵

        // 构造函数，初始化各个矩阵
        public PiResult(double[][] Pi, double[][] dPi, double[][] d2Pi, double[][] d3Pi) {
            this.Pi = Pi;
            this.dPi = dPi;
            this.d2Pi = d2Pi;
            this.d3Pi = d3Pi;
        }
    }

    // 计算项目反应理论 (IRT) 中的 Pi 矩阵及其导数
    public static PiResult Pi(double th, double[][] it, double D) {
        // 将项目参数矩阵转换为 RealMatrix 对象
        RealMatrix itemParams = new Array2DRowRealMatrix(it);
        // 获取项目的数量
        int nItems = itemParams.getRowDimension();
        // 获取每个项目的最大类别数（减去区分度参数）
        int maxCategories = itemParams.getColumnDimension();

        // 初始化概率矩阵及其导数矩阵
        double[][] prov = new double[nItems][maxCategories + 1];
        double[][] prov1 = new double[nItems][maxCategories + 1];
        double[][] prov2 = new double[nItems][maxCategories + 1];
        double[][] prov3 = new double[nItems][maxCategories + 1];

        // 遍历每个项目
        for (int i = 0; i < nItems; i++) {
            // 获取区分度参数 aj
            double aj = itemParams.getEntry(i, 0);
            // System.out.println("aj: " + aj);
            // 获取难度参数 bj
            double[] bj = itemParams.getSubMatrix(i, i, 1, maxCategories - 1).getData()[0];
            // System.out.println("bj: " + Arrays.toString(bj));

            // 初始化 Pjs 数组，包括边界条件
            double[] Pjs = new double[bj.length + 2];
            Pjs[0] = 1.0; // 边界条件：0分及0分以上的概率为1

            // 计算每个类别的概率 Pj
            for (int j = 0; j < bj.length; j++) {
                // System.out.println("bj: " + Arrays.toString(bj));
                double ej = Math.exp(D * aj * (th - bj[j]));
                Pjs[j +1] = ej / (1 + ej); // 根据公式计算中间概率
                // System.out.println("Pjs: " + Arrays.toString(Pjs));
            }

            Pjs[Pjs.length - 1] = 0.0; // 边界条件：大于最大分值的概率为0

            // System.out.println("Pjs: " + Arrays.toString(Pjs));


            // 初始化导数数组
            double[] dPjs = new double[Pjs.length];
            double[] d2Pjs = new double[Pjs.length];
            double[] d3Pjs = new double[Pjs.length];

            // 计算每个类别的导数
            for (int j = 0; j < Pjs.length; j++) {
                dPjs[j] = D * aj * Pjs[j] * (1 - Pjs[j]);
                d2Pjs[j] = D * aj * (dPjs[j] - 2 * Pjs[j] * dPjs[j]);
                d3Pjs[j] = D * aj * (d2Pjs[j] - 2 * Math.pow(dPjs[j], 2) - 2 * Pjs[j] * d2Pjs[j]);
            }

            // 计算相邻类别之间的差值
            for (int j = 0; j < maxCategories; j++) {
                prov[i][j] = Pjs[j] - Pjs[j + 1];
                prov1[i][j] = dPjs[j] - dPjs[j + 1];
                prov2[i][j] = d2Pjs[j] - d2Pjs[j + 1];
                prov3[i][j] = d3Pjs[j] - d3Pjs[j + 1];
            }
            // System.out.println("prov: " + Arrays.toString(prov));
        }

        // 返回计算结果
        return new PiResult(prov, prov1, prov2, prov3);
    }

    // 使用梯形积分法进行数值积分
    private static double intergrateCat(double[] x, double[] y) {
        // 创建梯形积分器
        TrapezoidIntegrator integrator = new TrapezoidIntegrator();
        // 定义被积函数
        UnivariateFunction function = t -> interpolate(x, y, t);
        // 执行积分
        // 它接受最大迭代次数、要积分的函数、积分区间的起点和终点作为参数，并返回积分的结果
        return integrator.integrate(2000, function, x[0], x[x.length - 1]);
    }

    // 线性插值函数
    private static double interpolate(double[] x, double[] y, double t) {
        // 处理边界情况
        if (t <= x[0]) return y[0];
        if (t >= x[x.length - 1]) return y[y.length - 1];
        // 进行线性插值
        for (int i = 0; i < x.length - 1; i++) {
            if (x[i] <= t && t <= x[i + 1]) {
                return y[i] + (y[i + 1] - y[i]) * (t - x[i]) / (x[i + 1] - x[i]);
            }
        }
        // 抛出异常，如果插值点超出范围
        throw new IllegalArgumentException("Interpolation point out of range.");
    }

    // 计算 EAP 估计值
    public static double eapEstGRM(double[][] it, int[] x, double D) {
        // 创建正态分布对象
        NormalDistribution normalDist = new NormalDistribution();

        // 定义对数似然函数 LL
        UnivariateFunction LL = th -> {
            PiResult piRes = Pi(th, it, D);
            double res = 1.0;
            // 计算每个项目的似然值
            for (int i = 0; i < x.length; i++) {
                res *= piRes.Pi[i][x[i]];
            }
            return res;
        };

        // 定义被积函数 g 和 h
        UnivariateFunction g = s -> s * normalDist.density(s) * LL.value(s);
        UnivariateFunction h = s -> normalDist.density(s) * LL.value(s);

        // 生成积分区间
        double[] X = linspace(-4, 4, 33);
        // 计算 g 和 h 在 X 上的值
        double[] Y1 = evaluate(X, g);
        double[] Y2 = evaluate(X, h);

        // 计算 EAP 估计值
        return intergrateCat(X, Y1) / intergrateCat(X, Y2);
    }

    // 生成等间距的数组
    private static double[] linspace(double start, double end, int num) {
        double[] result = new double[num];
        double step = (end - start) / (num - 1);
        for (int i = 0; i < num; i++) {
            result[i] = start + i * step;
        }
        return result;
    }

    // 计算函数在给定点上的值
    private static double[] evaluate(double[] x, UnivariateFunction f) {
        double[] result = new double[x.length];
        for (int i = 0; i < x.length; i++) {
            result[i] = f.value(x[i]);
        }
        return result;
    }


}
